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Modelling  of  Plane  Strain  Interfacial  Fracture  in  Incompressible  Materials 
T.  C.  Miller 

Air  Force  Research  Laboratory,  10  East  Saturn  Boulevard, 

Edwards  Air  Force  Base,  California  93524 
(805)  275-5323,  FAX  (805)  275-5435 

ABSTRACT 

Numerical  modelling  of  a  photoelastic  experiment  is  discussed.  The  experiment 
examined  incompressible  materials  under  plane  strain  conditions,  which  results  in  a  simplified 
analysis  due  to  a  vanishing  of  the  bimaterial  parameter.  The  photoelastic  experiment  used  the 
stress  freezing  method  to  determine  near  tip  stresses  in  interfacial  cracks  in  bimaterial  specimens. 
Different  crack  orientations  were  used  to  produce  different  mode  mixities.  Photoelastic  fringe 
patterns  were  analyzed  to  determine  the  magnitude  and  phase  angle  of  the  complex  stress 
intensity  factor.  These  experiments  were  modeled  using  a  finite  element  analysis  to  determine 
the  field  variables  near  the  tips  of  the  interfacial  cracks.  Magnitudes  of  the  complex  stress 
intensity  factors  are  found  from  J  integral  values,  derived  using  the  domain  integral  approach, 
and  the  phase  angles  are  determined  using  extrapolation  of  the  bond  line  traction  data  to_r  =  0. 
The  results  show  that  this  approach  is  a  useful  way  to  characterize  completely  the  complex  stress 
intensity  factor  in  incompressible  linear  elastic  bimaterial  combinations  under  plane  strain 
conditions. 

(Keywords:  interface,  fracture,  defects) 

INTRODUCTION 

The  solution  of  the  general  case  of  an  interfacial  crack  lying  between  two  linear  elastic 
isotropic  materials  was  first  introduced  by  Williams  in  1959’.  This  boundary  value  problem  was 
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solved  by  determining  the  eigenvalues  for  the  characteristic  equation,  resulting  in  an  infinite 
series  of  eigenvalues  with  imaginary  parts  related  to  the  material  properties  of  the  two  materials. 
The  imaginary  component,  €,  is  an  inherent  part  of  the  field  equations,  and  causes  a  number  of 
controversies  relating  to  inherent  mode  mixity,  stress  oscillation,  and  predictions  of  crack  face 
interpenetration2,3.  Later  works  attempted  to  explain  these  discrepancies4,5;  however,  in  a  large 
class  of  problems  these  contradictions  are  predicted  only  in  regions  very  close  to  the  crack  tip 
and  are  not  significant2.  Other  works  presented  after  Williams’  showed  that  other  eigenvalues 
existed  and  solved  the  problem  using  alternate  methods,  such  as  the  complex  variables 
formulation  of  Muskhelishvili.  Greater  detail  is  provided  in  several  sources2,3. 

A  subset  of  the  general  interfacial  fracture  problem  is  one  in  which  both  materials  are 
incompressible  and  plane  strain  conditions  exist.  In  this  case,  the  imaginary  part  of  the  complex 
eigenvalues  vanishes,  and  the  mechanics  of  the  crack  are  analogous  to  those  of  a  mixed-mode 
crack  in  a  homogenous  material6,7.  This  set  of  problems  includes  that  of  an  interfacial  crack  lying 
between  solid  rocket  propellant  and  a  rubber  liner.  To  investigate  the  propellant-liner 
relationship,  a  set  of  photoelastic  experiments  was  performed  that  also  used  incompressible 
materials6,8. 

In  previous  developments,  the  stress  intensity  factor  in  a  photoelastic  material  was  found 
using  the  set  of  isochromatic  fringe  loops  near  the  tip  of  the  crack.  The  location  of  the  point 
along  each  fringe  that  was  farthest  from  the  crack  tip  was  determined;  this  set  of  coordinate  data 
was  used  to  determine  the  stress  intensity  factor  by  extrapolating  the  data  to  r_=  0.  This 
technique  was  developed  for  mode  I  and  for  mixed-mode  loading9;  the  latter  was  used  for  the 
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photoelastic  experiments  discussed  here6.  Because  of  the  zero  valued  bimaterial  parameter^  the 
presence  of  two  materials  provides  no  additional  complications  when  using  this  procedure. 

The  photoelastic  experiments  used  single  edge  notched  tension  specimens  that 
incorporated  different  crack  orientations  to  produce  different  mode  mixities.  For  these  types  of 
specimen  and  boundary  conditions,  an  approximate  analytical  solution  could  be  constructed6. 
However,  more  complex  geometries  preclude  the  use  of  an  analytical  expression  for  stress 
intensity  factor  calculations.  Also,  many  bimaterial  combinations  cannot  be  analyzed  using 
analogous  photoelastic  materials  without  ignoring  important  features  of  the  material  behavior.  In 
such  cases,  developing  numerical  models  may  confirm  experimental  results  and  help  extend  the 
predictive  capabilities  of  the  research.  It  is  with  this  goal  that  the  numerical  models  featured  here 
have  been  developed.  These  models  simulate  the  photoelastic  tests  so  that  the  evaluation  of 
stress  intensity  factor  magnitudes  and  phase  angles  from  numerical  models  can  be  assessed.  A 
similar  method  will  be  used  in  the  study  of  incompressible  bimaterial  combinations  related  to 
solid  rocket  motor  design. 

THEORY 

In  general,  an  interfacial  crack  between  two  isotropic  linear  elastic  materials  will  be 
characterized  by  a  bimaterial  parameter  erk 
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Here  P  is  one  of  the  two  Dundurs'  elastic  mismatch  parameters,  p  is  the  shear  modulus, _v 
is  Poisson’s  ratio,  and  k  =  (3  -  v)  /  (1  +  v)  in  plane  stress  and  (3  -  4v)  in  plane  strain.  The 
subscripts  are  used  to  index  the  two  materials.  In  a  region  near  the  crack  tip,  the  higher  order 
terms  are  negligible  and  the  field  expressions  depend  only  on  the  singular  terms10: 

a  ‘  —{Re(Kr“i£„(ff)  *  Im(Kr“)^q<0)}  (2) 

\[2nr 

Here  the  two  S  functions  characterize  the  angular  variations  of  the  near  tip  field. 

<v\ 

Because  in-plane  tensile  and  shear  loading  are  coupled,  the  parameter  K  characterizes  the 
combination  of  these  two  effects  as  a  complex  stress  intensity  factor.  The  inherent  mode  mixity 
is  manifested  in  the  expression  for  the  traction  along  the  interface,  derived  from  eqn  (2)  by  letting 
£=  02’3: 

( o  +iojg-0  =  [cos(eLnr)  +  isin(eLnr)]  (3) 

^  ^  sfTrcr  JJTr 

The  complex  stress  intensity  factorK  =  K,  +  i  K2  =  Ke1'1'  depends  on  the  geometry,  loads, 
stress  state,  and  materials.  Inspection  of  eqn  (3)  shows  that  the  bond  line  traction  ratio 
(oxy  /  o>T)e = 0  varies  with  distance  from  the  crack  tip.  In  a  region  near  the  crack  tip,  the  phase 
angles  of  the  left-hand  and  right-hand  sides  of  eqn  (3)  differ  by  e  Ln  r: 

tan~'[(ax/ayy)d=(J  =  tan'1  [K/K,]  +  eLnr  (4) 

The  magnitude  of  K  is  related  to  J,  the  contour  integral,  through2: 

j  _  (c,  *  cj  [Kf_ 

16  cosh2(ne) 


(5) 
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Here  cp  =  4(1  -  vp)  /  pp  in  plane  strain  and  4  /  (pp(l  +  vp))  in  plane  stress.  Equations  (1-5) 
show  that  the  presence  of  the  nonzero  bimaterial  parameter^  causes  complexities  not  present  in 
the  consideration  of  homogeneous  materials.  However,  in  the  special  case  where  both  materials 
are  incompressible  and  plane  strain  conditions  prevail,  then  v,  =  v2  =  14  and  Kp  =  (3  -  4vp); 
substitution  into  eqn  (1)  shows  that J3_=  _e_=  0  in  these  circumstances.  Introducing  e_=  0  into 
eqns  (3-5)  gives  simplified  expressions  for  the  degenerate  case6,7:* 


(6) 


The  phase  angle  of  the  complex  stress  intensity  factor  is  shown  here  as  Y  and  the 
effective  plane  strain  elastic  modulus  E*  depends  on  the  plane  strain  elastic  moduli  of  the  two 
materials.  Equations  (6-8)  show  that  the  solution  for  an  interfacial  crack  in  an  incompressible 
bimaterial  pair  under  plane  strain  conditions  is  analogous  to  that  of  a  homogeneous  material  with 

*  Because  v,  =  v2  =  1/2,  additional  simplifications  could  be  made  to  eqn  (8),  but  it  has 
been  presented  as  shown  to  preserve  the  meaning  of  the  parameters. 
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mixed-mode  loading.  The  ratios  (Kj  /  K,)  and  (oxy  /  o^.  )e=0  are  equal  in  the  near  tip  region,  andj_ 
and  K  are  related  through  a  modulus  parameter.  The  additional  complexities  of  stress 
oscillations,  inherent  mode  mixity,  and  crack  face  interpenetration  are  no  longer  present. 

In  the  photoelastic  experiment,  mode  mixity  was  varied  by  maintaining  the  vertical 
loading  direction  while  testing  specimens  that  have  different  crack  orientations.  Denoting  the 
angle  the  crack  makes  with  respect  to  the  mode  I  loading  orientation  by  T,  the  photoelastic 
experiments  used  specimens  with  crack  orientations  of T  =  0°,  15°,  30°,  and  45°.  ForT  =  0° 
(mode  I  loading),  the  stress  intensity  factor  phase  angle  equals  zero.  The  mode  mixity  and  phase 
angle  of  K  increase  in  an  approximately  linear  way  with  increases  in  F5'8. 

DISCUSSION 

Figure  1  shows  a  representative  specimen  modelled  using  finite  element  analysis6  8.  The 
set  of  investigations  on  which  the  numerical  models  is  based  used  a  photoelastic  polymer, 
araldite,  loaded  above  its  stress  freezing  temperature  of  1 16°  C.  The  second  material  is 
composed  of  araldite  and  aluminum  powder,  so  that  the  two  materials  have  moduli  of  18.6  and 
36.9  MPa,  respectively,  above  the  stress  freezing  temperature.  The  latter  material  is  opaque  due 
to  the  introduction  of  the  aluminum  powder,  so  that  photoelastic  fringe  data  is  only  available  for 
the  upper  half  of  each  specimen.  Above  the  stress  freezing  temperature,  both  materials  are 
incompressible.  The  two  materials  were  bonded  together  using  a  thin  layer  (less  than  0.5  mm)  of 
photoelastic  adhesive.  The  material  properties  of  the  adhesive  layer  were  similar  to  that  of  the 
araldite-aluminum;  later  work  confirmed  that  the  presence  of  a  thin  adhesive  layer  causes  no 
significant  changes  in  the  relationship  between  the  loads  and  the  complex  stress  intensity  factor  . 
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To  simulate  cracks,  notches  were  machined  along  the  interfaces  to  a  depth  of  9.5  mm. 
Loading  was  accomplished  using  freely  rotation  aluminum  grips  to  which  the  specimens  were 
bonded,  and  photoelastic  fringe  patterns  of  the  loaded  specimens  were  recorded  for  analysis7,8. 

The  complex  stress  intensity  factor  was  analyzed  using  fringe  loop  data  from  the  araldite 
portion  of  the  specimens.  The  fringe  orders  and  locations  of  the  maximum  radii  of  the  fringe 
loops  near  the  crack  tip  were  recorded.  By  extrapolating  the  data  to  the  crack  tip,  the  complex 
stress  intensity  factor  can  be  evaluated.  This  method  overestimated  the  value  of  applied  K  for  the 
smaller  crack  angles  because  of  the  influence  of  residual  stresses  and  machining  stresses  on  the 
fringe  patterns.  An  analytical  model  was  developed  that  modified  predictions  for  an  edge 
cracked  geometry  with  finite  width  effects  and  notch  effects.  The  experimental  results  shown  in 
this  work  are  the  results  of  analysis  of  stress  frozen  slices  taken  from  the  mid-thickness  of  the 
specimens  so  that  the  data  from  the  experiment  represents  plane  strain  conditions6,8. 

Figure  1  shows  a  typical  finite  element  mesh.  Eight-noded  quadrilateral  elements  are 
used  throughout  the  mesh.  Quarter  point  elements  surround  the  crack  tip  in  a  spider  web 
formation  as  shown  in  Figure  2.  The  boundary  conditions  were  vertical  displacements  applied  at 
the  outside  middle  nodes  of  the  aluminum  grips  (see  Figure  1 ).  This  type  of  loading  closely 
resembled  the  pin  loading  of  the  actual  specimens;  rotation  at  the  grips  significantly  affects  the 
mode  mixity  at  the  crack  tip. 

One  complication  that  arises  in  the  modelling  of  elastic  incompressible  materials  is  the 
indeterminacy  of  conventional  finite  element  formulations.  Typically,  a  finite  element 
formulation  solves  for  the  displacement  components  at  the  nodes  in  the  mesh.  The  solution  will 
be  the  configuration  that  minimizes  the  potential  energy  of  the  mesh,  and,  in  the  approximate 
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sense,  the  potential  energy  of  the  continuum  that  the  mesh  represents.  The  strains  and  stresses 
are  derived  from  the  nodal  displacements  using  shape  functions  and  constitutive  relationships. 
Usually,  a  single  set  of  displacements,  strains,  and  stresses  minimizes  the  potential  energy  of  the 
mesh.  However,  with  incompressible  materials,  addition  of  an  arbitrary  pressure  term  yields 
another  solution  that  minimizes  the  potential  energy,  so  that  there  are  an  infinite  number  of 
solutions.  To  resolve  this  dilemma,  the  finite  element  problem  is  formulated  with  the  hydrostatic 
component  of  stress  as  an  additional  solution  variable.  For  the  eight-noded  quadrilateral 
elements  used  here,  this  resulted  in  three  additional  degrees  of  freedom  per  element:  a  constant 
term  and  linear  variations  in  two  orthogonal  directions.  The  resulting  mixed  formulation  avoids 
the  static  indeterminacy  associated  with  purely  displacement-based  formulations12'14. 

The  results  for  each  specimen  include  an  estimate  of  the 2  integral  based  on  a  domain 
integral  approach.  To  implement  this  approach,  two  concentric  contours  are  described  that  begin 
on  the  lower  crack  face  and  end  on  the  upper  one.  The_J  integral  for  the  inner  contour  may  be 
restated  as  an  integral  involving  both  contours  by  using  a  smooth  weighting  function15,16: 

J  =  / (au-^r  "  w5^] mjq i dc  (10) 

c  1 

Here  the  closed  contour  C  is  composed  of  four  contours:  C,nner,  the  inner  contour;  C01iter, 
the  outer  contour,  Cupper,  the  portion  of  the  upper  crack  face  (see  Figure  3),  and  C,ower,  the 
corresponding  portion  of  the  lower  crack  face.  Also,  w  is  the  strain  energy  density,  nij  is  the  jth 
component  of  the  normal  pointing  away  from  the  enclosed  area,  and  q,  is  the  x,  component  of  a 
certain  smoothing  function,  q.  The  above  equation  can  be  converted  to  an  equivalent  area 
integral  using  the  Gauss  divergence  theorem,  giving: 


9 


J  = 


du 
(a.. — - 
,Jdx. 

A  1 


/° 


<ii) 


Here  A_  is  the  area  enclosed  by  the  contours.  In  practice,  the  inner  contour  is  allowed  to 
shrink  onto  the  crack  tip,  the  outer  contour  extends  along  the  edges  of  the  quadrilateral  elements, 
and  the  area  integral  is  approximated  using  Gauss  quadrature  formulas  and  values  for  the  field 
variables  evaluated  at  the  integration  points  for  each  element16  '7.  A  method  for  evaluating  J  for 
general  interfacial  cracks  in  bimaterials  is  given  by  Shih  and  Asaro,  and  allow  for  separation  of 
the  energy  contributions  and  determination  of  K,  and  K2  using  an  interaction  energy  release 
rate16.  For  the  degenerate  case  discussed  here,  a  simpler  methodology  can  be  used  to  determine 
the  overall  energy  release  rate15.  The  domain  integral  method  of  estimating  J  works  well  with 
finite  element  methods  and  provides  for  a  robust  determination  of  J  that  is  not  sensitive  to  mesh 
construction  and  that  does  not  require  highly  accurate  stress  approximations  very  near  to  the 
crack  tip.  Also,  many  finite  element  software  programs  have  built-in  algorithms  for  the 
determination  of  the  J_  integral  in  this  manner.  The  J  integral  value  was  used  with  eqn  (8)  to 
calculate  the  magnitude,  K,  of  the  complex  stress  intensity  factor. 

Accurate  evaluation  of  the  phase  angle,  *P,  is  then  required.  Equations  (6)  and  (7)  show 
that  the  arctangent  of  the  traction  ratio  ojo^  along  the  bond  line  is  equal  to  Y  in  the  near  tip 
region.  Using  curve  fitting  techniques,  tan'1[(oxy/ayy)9=0  ]  can  be  evaluated  as  r  -  0  and  taken  as 
an  approximation  of Y  (for  the  work  performed  here,  a  cubic  polynomial  was  used).  One 
problem  with  this  method  is  that  inaccuracies  in  the  stress  components  caused  by  high  gradients . 
near  the  crack  tip  make  extrapolation  of  (oxy/oyy)e =0  to_r  =  0  difficult.  An  unrefined  mesh  that 
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incorporates  only  a  few  quarter  point  elements  will  exhibit  excessive  scatter  in  tan'1[(oxy/oyy)0=o  ] 
as_r  -  0,  so  that  this  regression  technique  is  unfeasible.  This  method  can  be  implemented, 
however,  with  a  more  refined  mesh,  such  as  the  one  used  here,  which  incorporates  approximately 
30  quarter  point  elements  (the  number  differs  for  the  various  crack  geometries  considered). 

The  stresses  in  a  finite  element  formulation  are  typically  found  using  derivatives  of 
displacement  fields  based  on  nodal  displacements.  Consequently,  results  for  stresses  are 
generally  less  accurate  than  displacements  in  any  finite  element  solution.  This  dictates  that, 
when  possible,  extrapolation  techniques  should  be  applied  to  displacements  rather  than  stresses. 
A  method  analogous  to  the  stress-based  extrapolation  technique  exists  for  relative  crack  face 
displacements1018.  However,  this  method  requires  that  the  two  crack  faces  be  initially  coincident, 
so  that  this  technique  cannot  be  easily  applied  to  the  current  geometry,  which  incorporates  crack 
faces  with  finite  separation. 

In  this  work,  J  integral  values  were  used  to  evaluate  K,  and  bond  line  traction  data  was 
used  to  determine  Y,  the  phase  angle  of  K.  Results  are  presented  and  compared  with  the 

-  A S\J 

photoelasticity  results  next. 

RESULTS 

Figure  4  shows  a  set  of  representative  contour  plots;  the  plots  are  for  a  crack  orientation 
of  45°.  The  relative  size  of  the  two  sets  of  fringe  loops  is  caused  by  mismatch  of  the  two  elastic 
moduli;  the  shape  of  the  o}y  and  oxy  contours  remains  invariant  with  respect  to  mode  mixity. 

Other  stress  variables,  however,  such  as  the  maximum  in-plane  stress,  are  strong  functions  of 
mode  mixity.  Table  1  summarizes  the  results  for  both  the  magnitude  and  phase  of  K. 
Experimental  and  computational  results  agreed  well.  With  the  T  =  0°  and  T  =15°  specimens,  the 
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photoelastic  results  were  high  due  to  the  effect  previously  mentioned:  residual  and  machining 
stresses  incorporated  into  the  specimen  during  fabrication  influence  the  gradient  of  the  fringe 
patterns  in  the  vertical  direction,  and  this  gradient  is  used  in  the  determination  of  K. 

Table  1  also  shows  the  phase  angle  results.  The  worst  discrepancy  between  the  numerical 
and  experimental  results  was  for  the  T  =  45°  bimaterial  specimen,  in  which  case  the  values 
differed  by  4.2°.  This  is  accurate  enough  for  most  applications,  however,  if  necessary,  other 
methods  may  be  used  to  obtain  additional  improvements  in  accuracy  of  phase  angle 
determination' 6’17.  However,  the  method  shown  above  is  often  preferable  because  of  its 
simplicity  compared  with  other  procedures. 

CONCLUSIONS 

The  fracture  mechanics  of  interfacial  cracks  in  incompressible  bimaterials  subject  to 
plane  strain  conditions  closely  resembles  those  for  homogeneous  materials  subject  to 
mixed-mode  loading.  Static  indeterminacy  can  cause  inaccuracies  unless  a  mixed  formulation  is 
used.  Numerical  modelling  can  be  used  to  characterize  the  complex  stress  intensity  factors  for 
various  specimen  geometries.  The  magnitude  ofK^can  be  evaluated  from  J_  integral  calculations, 
and  the  phase  angle  of  K  can  be  determined  by  extrapolating  field  variables  such  as  the  bond  line 
traction  data.  Mesh  refinements  should  include  numerous  quarter  point  elements  near  the  crack 
tip  so  that  phase  angles  can  be  evaluated  using  extrapolation  of  bond  line  traction  data. 
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Table  1  Comparison  of  numerical  modelling  and  photoelastic 
results  for  stress  intensity  factors 
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computational 
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Abstract 

Numerical  modelling  of  a  photoelastic  experiment  is  discussed.  The  experiment  examined  incompressible  materials  under  plane  strain 
conditions,  which  results  in  a  simplified  analysis  owing  to  a  vanishing  of  the  bimaterial  parameter.  The  photoelastic  experiment  used  the 
stress  freezing  method  to  determine  near  tip  stresses  in  interfacial  cracks  in  bimaterial  specimens.  Different  crack  orientations  were  used  to 
produce  different  mode  mixities.  Photoelastic  fringe  patterns  were  analyzed  to  determine  the  magnitude  and  phase  angle  of  the  complex 
stress  intensity  factor.  These  experiments  were  modeled  using  a  finite  element  analysis  to  determine  the  field  variables  near  the  tips  of  the 
interfacial  cracks.  Magnitudes  of  the  complex  stress  intensity  factors  are  found  from  J  integral  values,  derived  using  the  domain  integral 
approach,  and  the  phase  angles  are  determined  using  extrapolation  of  the  bond  line  traction  data  to  r  =  0.  The  results  show  that  this  approach 
is  a  useful  way  to  characterize  completely  the  complex  stress  intensity  factor  in  incompressible  linear  elastic  bimaterial  combinations  under 
plane  strain  conditions.  ©  1999  Elsevier  Science  Ltd.  All  rights  reserved. 

Keywords:  A.  Interface/interphase;  B.  Fracture;  B.  Defects 


Nomenclature 

/3  Dundur’s  second  elastic  mismatch  parameter 

T  Crack  orientation  angle 

e  Bimaterial  parameter 

6  Angular  orientation  (second  polar  coordinate) 
k  Kappa 

fi  Shear  modulus 

v  Poisson’s  ratio 

t Tij  Stress  tensor  component 

'VP  Complex  stress  intensity  factor  phase  angle 

£*  Bimaterial  effective  plane  strain  modulus 

Ei,  E2  Plane  strain  moduli  for  component  materials 
J  J  integral 

K  Complex  stress  intensity  factor 

K  complex  stress  intensity  factor  magnitude 

r  radius  (first  polar  coordinate) 


1.  Introduction 

The  solution  of  the  general  case  of  an  interfacial  crack 
lying  between  two  linear  elastic  isotropic  materials  was  first 
introduced  by  Williams  in  1959  [1].  This  boundary  value 
problem  was  solved  by  determining  the  eigenvalues  for  the 
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characteristic  equation,  resulting  in  an  infinite  series  of 
eigenvalues  with  imaginary  parts  related  to  the  material 
properties  of  the  two  materials.  The  imaginary  component, 
e,  is  an  inherent  part  of  the  field  equations,  and  causes  a 
number  of  controversies  relating  to  inherent  mode  mixity, 
stress  oscillation,  and  predictions  of  crack  face  interpenetra¬ 
tion  [2,3].  Later  works  attempted  to  explain  these  discrepan¬ 
cies  [4,5];  however,  in  a  large  class  of  problems  these 
contradictions  are  predicted  only  in  regions  very  close  to 
the  crack  tip  and  are  not  significant  [2].  Other  works 
presented  after  Williams’  showed  that  other  eigenvalues 
existed  and  solved  the  problem  using  alternate  methods, 
such  as  the  complex  variables  formulation  of  Muskhelish- 
vili.  Greater  detail  is  provided  in  several  sources  [2,3]. 

A  subset  of  the  general  interfacial  fracture  problem  is  one 
in  which  both  materials  are  incompressible  and  plane  strain 
conditions  exist.  In  this  case,  the  imaginary  part  of  the 
complex  eigenvalues  vanishes,  and  the  mechanics  of  the 
crack  are  analogous  to  those  of  a  mixed-mode  crack  in  a 
homogenous  material  [6,7].  This  set  of  problems  includes 
that  of  an  interfacial  crack  lying  between  solid  rocket 
propellant  and  a  rubber  liner.  To  investigate  the  propel¬ 
lant-liner  relationship,  a  set  of  photoelastic  experiments 
was  performed  that  also  used  incompressible  materials 
[6,8]. 

In  previous  developments,  the  stress  intensity  factor  in  a 
photoelastic  material  was  found  using  the  set  of  isochro- 
matic  fringe  loops  near  the  tip  of  the  crack.  The  location 
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of  the  point  along  each  fringe  that  was  farthest  from  the 
crack  tip  was  determined;  this  set  of  coordinate  data  was 
used  to  determine  the  stress  intensity  factor  by  extrapolating 
the  data  to  r  =  0.  This  technique  was  developed  for  mode  I 
and  for  mixed-mode  loading  [9];  the  latter  was  used  for  the 
photoelastic  experiments  discussed  here  [6].  Because  of  the 
zero  valued  bimaterial  parameter  e,  the  presence  of  two 
materials  provides  no  additional  complications  when  using 
this  procedure. 

The  photoelastic  experiments  used  single  edge  notched 
tension  specimens  that  incorporated  different  crack  orienta¬ 
tions  to  produce  different  mode  mixities.  For  these  types  of 
specimen  and  boundary  conditions,  an  approximate  analy¬ 
tical  solution  could  be  constructed  [6].  However,  more 
complex  geometries  preclude  the  use  of  an  analytical 
expression  for  stress  intensity  factor  calculations.  Also, 
many  bimaterial  combinations  cannot  be  analyzed  using 
analogous  photoelastic  materials  without  ignoring  impor¬ 
tant  features  of  the  material  behavior.  In  such  cases,  devel¬ 
oping  numerical  models  may  confirm  experimental  results 
and  help  extend  the  predictive  capabilities  of  the  research.  It 
is  with  this  goal  that  the  numerical  models  featured  here 
were  developed.  These  models  simulate  the  photoelastic 
tests  so  that  the  evaluation  of  stress  intensity  factor  magni¬ 
tudes  and  phase  angles  from  numerical  models  can  be 
assessed.  A  similar  method  will  be  used  in  the  study  of 
incompressible  bimaterial  combinations  related  to  solid 
rocket  motor  design. 


2.  Theory 


In  general,  an  interfacial  crack  between  two  isotropic 
linear  elastic  materials  will  be  characterized  by  a  bimaterial 
parameter  e  [2,3]. 


1  inf  1  ~  p-  ~  ~  /A2(K|  ~ 

2tT  V  1  +  (3  /  jJL j(k2  +  1)  +  £t2(Kl  +  1) 


(1) 


Here  (3  is  one  of  the  two  Dundurs’  elastic  mismatch  para¬ 
meters,  fi  the  shear  modulus,  v  Poisson’s  ratio,  and  k  = 
(3  —  v)/(l  +  v)  in  plane  stress  and  (3  —  Av)  in  plane  strain. 
The  subscripts  are  used  to  index  the  two  materials.  In  a 
region  near  the  crack  tip,  the  higher  order  terms  are  negli¬ 
gible  and  the  field  expressions  depend  only  on  the  singular 
terms  [10]: 


V2 m 


1  U 

Re(tfri£)  £  (0)  +  £  (#) 

pq  pq 


(2) 


Here  the  two  X  functions  characterize  the  angular  variations 
of  the  near  tip  field.  Because  in-plane  tensile  and  shear 
loading  are  coupled,  the  parameter  K  characterizes  the 
combination  of  these  two  effects  as  a  complex  stress  inten¬ 
sity  factor.  The  inherent  mode  mixity  is  manifested  in  the 
expression  for  the  traction  along  the  interface,  derived  from 


Eq.  (2)  by  letting  (9=0  [2,3]: 

Krl€  K 


(crvv  +  io-vv)  =  —==  =  -=r[cos(eln  r)  +  /sin(elnr)] 

\  --  ■/((--{)  ^ 2m‘  yJATTP 

(3) 


The  complex  stress  intensity  factor  K  —  K{  +  iK2  = 

Ke1^  depends  on  the  geometry,  loads,  stress  state,  and  mate¬ 
rials.  Inspection  of  Eq.  (3)  shows  that  the  bond  line  traction 
ratio  (rrvv/rrvv)^0  varies  with  distance  from  the  crack  tip.  In 
a  region  near  the  crack  tip,  the  phase  angles  of  the  left-hand 
and  right-hand  sides  of  Eq.  (3)  differ  by  e  In  r : 

'anifeLh"1f]+e,n'  (4) 

The  magnitude  of  K  is  related  to  7,  the  contour  integral, 
through  [2]: 

j  =  <£i±s>II:  (5) 

16  COS/?2(7T6) 

P 

Here  cp  —  4(1  -  vp)!i±p  in  plane  strain  and  4/(/x/?(l  +  vp)) 
in  plane  stress.  Eqs.  (l)-(5)  show  that  the  presence  of  the 
nonzero  bimaterial  parameter  e  causes  complexities  not  i 
present  in  the  consideration  of  homogeneous  materials. 
However,  in  the  special  case  where  both  materials  are 
incompressible  and  plane  strain  conditions  prevail,  then 
v\  =  v2  =  1/2  and  Kp  —  (3  -  Avp)\  substitution  into  Eq. 

(1)  shows  that  [3  =  e  —  0  in  these  circumstances.  Introdu¬ 
cing  €  =  0  into  Eqs.  (3)~(5)  gives  simplified  expressions  for 


the  degenerate  case  [6,7]  : 

K 

((7VV  +  l  crvv.)^0  —  , 

V  2777’ 

(6) 

)  ] 

AV  /  0=0 J 

(7) 

E* 

E*  2  L  £,  E2 

The  phase  angle  of  the  complex  stress  intensity  factor  is 
shown  here  as  and  the  effective  plane  strain  elastic  modu¬ 
lus  E*  depends  on  the  plane  strain  elastic  moduli  of  the  two 
materials.  Eqs.  (6)-(8)  show  that  the  solution  for  an  inter¬ 
facial  crack  in  an  incompressible  bimaterial  pair  under  plane 
strain  conditions  is  analogous  to  that  of  a  homogeneous 
material  with  mixed-mode  loading.  The  ratios  (K2/K{)  and 
(crvv/crvv)0=o  are  equal  in  the  near  tip  region  and  J  and  K  are 
related  through  a  modulus  parameter.  The  additional 


£i  = 


(8) 


i-*E 


E2  = 


A 


1  Because  v\  —  v2=  1/2,  additional  simplifications  could  he  made  to  Eq. 
(8),  but  it  has  been  presented  as  shown  to  preserve  the  meaning  of  the 
parameters. 
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Fig.  1.  Typical  finite  element  mesh  used  in  numerical  computations  (crack 
orientation  =15  degrees). 

complexities  of  stress  oscillations,  inherent  mode  mixity, 
and  crack  face  interpenetration  are  no  longer  present. 

In  the  photoelastic  experiment,  mode  mixity  was  varied 
by  maintaining  the  vertical  loading  direction  while  testing 
specimens  that  have  different  crack  orientations.  Denoting 
the  angle  the  crack  makes  with  respect  to  the  mode  I  loading 
orientation  by  T,  the  photoelastic  experiments  used 


Fig.  2.  Detail  of  finite  element  mesh  near  the  crack  tip  (crack  orientation  — 
15  degrees). 


specimens  with  crack  orientations  of  T  =  0°,  15°,  30°,  and 
45°.  For  T  =  0°  (model  I  loading),  the  stress  intensity  factor 
phase  angle  equals  zero.  The  mode  mixity  and  phase  angle 
of  K  increase  in  an  approximately  linear  way  with  increases 
in  T  [6,8]. 


3.  Discussion 

Fig.  1  shows  a  representative  specimen  modelled  using 
finite  element  analysis  [6,8].  The  set  of  investigations  on 
which  the  numerical  models  is  based  used  a  photoelastic 
polymer,  araldite,  loaded  above  its  stress  freezing  tempera¬ 
ture  of  1 16°C.  The  second  material  is  composed  of  araldite 
and  aluminum  powder,  so  that  the  two  materials  have 
moduli  of  18.6  and  36.9  MPa,  respectively,  above  the  stress 
freezing  temperature.  The  latter  material  is  opaque  because 
of  the  introduction  of  the  aluminum  powder,  so  that  photo¬ 
elastic  fringe  data  is  only  available  for  the  upper  half  of  each 
specimen.  Above  the  stress  freezing  temperature,  both 
materials  are  incompressible.  The  two  materials  were 
bonded  together  using  a  thin  layer  (less  than  0.5  mm)  of 
photoelastic  adhesive.  The  material  properties  of  the  adhe¬ 
sive  layer  were  similar  to  that  of  the  araldite- aluminum; 
later  work  confirmed  that  the  presence  of  a  thin  adhesive 
layer  causes  no  significant  changes  in  the  relationship 
between  the  loads  and  the  complex  stress  intensity  factor 
[11]. 

To  simulate  cracks,  notches  were  machined  along  the 
interfaces  to  a  depth  of  9.5  mm.  Loading  was  accomplished 
using  freely  rotating  aluminum  grips  to  which  the  specimens 
were  bonded,  and  photoelastic  fringe  patterns  of  the  loaded 
specimens  were  recorded  for  analysis  [7,8]. 

The  complex  stress  intensity  factor  was  analyzed  using 
fringe  loop  data  from  the  araldite  portion  of  the  specimens. 
The  fringe  orders  and  locations  of  the  maximum  radii  of  the 
fringe  loops  near  the  crack  tip  were  recorded.  By  extrapo¬ 
lating  the  data  to  the  crack  tip,  the  complex  stress  intensity 
factor  can  be  evaluated.  This  method  overestimated  the 
value  of  applied  K  for  the  smaller  crack  angles  because  of 
the  influence  of  residual  stresses  and  machining  stresses  on 
the  fringe  patterns.  An  analytical  model  was  developed  that 
modified  predictions  for  an  edge  cracked  geometry  with 
finite  width  effects  and  notch  effects.  The  experimental 
results  shown  in  this  work  are  the  results  of  analysis  of  stress 
frozen  slices  taken  from  the  mid-thickness  of  the  specimens 
so  that  the  data  from  the  experiment  represents  plane  strain 
conditions  [6,8]. 

Fig.  1  shows  a  typical  finite  element  mesh.  Eight-noded 
quadrilateral  elements  are  used  throughout  the  mesh.  Quar¬ 
ter  point  elements  surround  the  crack  tip  in  a  spider  web 
formation  as  shown  in  Fig.  2.  The  boundary  conditions  were 
vertical  displacements  applied  at  the  outside  middle  nodes 
of  the  aluminum  grips  (see  Fig.  1).  This  type  of  loading 
closely  resembled  the  pin  loading  of  the  actual  specimens; 
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y  material  1 


material  2 


Fig.  3.  Evaluation  of  J  integral  using  domain  integral  method. 

rotation  at  the  grips  significantly  affects  the  mode  mixity  at 
the  crack  tip. 

One  complication  that  arises  in  the  modelling  of  elastic 
incompressible  materials  is  the  indeterminacy  of  conven¬ 
tional  finite  element  formulations.  Typically,  a  finite 
element  formulation  solves  for  the  displacement  compo¬ 
nents  at  the  nodes  in  the  mesh.  The  solution  will  be  the 
configuration  that  minimizes  the  potential  energy  of  the 
mesh,  and,  in  the  approximate  sense,  the  potential  energy 
of  the  continuum  that  the  mesh  represents.  The  strains  and 
stresses  are  derived  from  the  nodal  displacements  using 
shape  functions  and  constitutive  relationships.  Usually,  a 
single  set  of  displacements,  strains,  and  stresses  minimizes 
the  potential  energy  of  the  mesh.  However,  with  incompres¬ 
sible  materials,  addition  of  an  arbitrary  pressure  term  yields 
another  solution  that  minimizes  the  potential  energy,  so  that 
there  are  an  infinite  number  of  solutions.  To  resolve  this 
dilemma,  the  finite  element  problem  is  formulated  with 
the  hydrostatic  component  of  stress  as  an  additional  solution 
variable.  For  the  eight-noded  quadrilateral  elements  used 
here,  this  resulted  in  three  additional  degrees  of  freedom 
pet  element:  a  constant  term  and  linear  variations  in  two 
orthogonal  directions.  The  resulting  mixed  formulation 
avoids  the  static  indeterminacy  associated  with  purely 
displacement-based  formulations  [12-14]. 

The  results  for  each  specimen  include  an  estimate  of  the  J 
integral  based  on  a  domain  integral  approach.  To  implement 
this  approach,  two  concentric  contours  are  described  that 
begin  on  the  lower  crack  face  and  end  on  the  upper  one. 
The  J  integral  for  the  inner  contour  may  be  restated  as  an 
integral  involving  both  contours  by  using  a  smooth  weight¬ 
ing  function  [15,16]. 

J  =  \c{ar‘jtoi  ~  w8u)mj<J  1  dc  (9) 

Here  the  closed  contour  C  is  composed  of  four  contours: 
Cinner,  the  inner  contour,  Coutei>  the  outer  contour,  Cupper,  the 
portion  of  the  upper  crack  face  (see  Fig.  3)  and  C!ovver,  the 
corresponding  portion  of  the  lower  crack  face.  Also,  w  is  the 
strain  energy  density,  m}  is  the  yth  component  of  the  normal 
pointing  away  from  the  enclosed  area,  and  <71  is  the  X\ 


component  of  a  certain  smoothing  function,  q .  The  above4 
equation  can  be  converted  to  an  equivalent  area  integral 
using  the  Gauss  divergence  theorem,  giving: 


Here  A  is  the  area  enclosed  by  the  contours.  In  practice,  the 
inner  contour  is  allowed  to  shrink  onto  the  crack  tip,  the 
outer  contour  extends  along  the  edges  of  the  quadrilateral 
elements,  and  the  area  integral  is  approximated  using  Gauss 
quadrature  formulas  and  values  for  the  field  variables  eval¬ 
uated  at  the  integration  points  for  each  element  [16,17].  A 
method  for  evaluating  J  for  general  interfacial  cracks  in 
bimaterials  is  given  by  Shih  and  Asaro,  and  allows  for 
separation  of  the  energy  contributions  and  determination 
of  K\  and  K2  using  an  interaction  energy  release  rate  [16]. 
For  the  degenerate  case  discussed  here,  a  simpler  methodol¬ 
ogy  can  be  used  to  determine  the  overall  energy  release  rate 
[15].  The  domain  integral  method  of  estimating  J  works 
well  with  finite  element  methods  and  provides  for  a  robust 
determination  of  J  that  is  not  sensitive  to  mesh  construction 
and  that  does  not  require  highly  accurate  stress  approxima¬ 
tions  very  near  to  the  crack  tip.  Also,  many  finite  element 
software  programs  have  built-in  algorithms  for  the  determi¬ 
nation  of  the  J  integral  in  this  manner.  The  J  integral  value 
was  used  with  Eq.  (8)  to  calculate  the  magnitude,  K ,  of  the 
complex  stress  intensity  factor. 

Accurate  evaluation  of  the  phase  angle,  3P,  is  then 
required.  Eqs.  (6)  and  (7)  show  that  the  arctangent  of  the 
traction  ratio  axyJayy  along  the  bond  line  is  equal  to  XP  in  the 
near  tip  region.  Using  curve  fitting  techniques, 
tan^^cr^/a^O^ol  can  be  evaluated  as  r  — ►  0  and  taken 
as  an  approximation  of  (for  the  work  performed  here,  a 
cubic  polynomial  was  used).  One  problem  with  this  method 
is  that  inaccuracies  in  the  stress  components  caused  by  high 
gradients  near  the  crack  tip  make  extrapolation  of 
(rrvv/crvv)0=o  to  r  =  0  difficult.  An  unrefined  mesh  that  incor¬ 
porates  only  a  few  quarter  point  elements  will  exhibit  exces¬ 
sive  scatter  in  tan-1[(crvv/ayv)f?=0]  as  /  *— ►  0,  so  that  this 
regression  technique  is  unfeasible.  This  method  can  be 
implemented,  however,  with  a  more  refined  mesh,  such  as 
the  one  used  here,  which  incorporates  approximately  30 
quarter  point  elements  (the  number  differs  for  the  various 
crack  geometries  considered). 

The  stresses  in  a  finite  element  formulation  are  typically 
found  using  derivatives  of  displacement  fields  based  on 
nodal  displacements.  Consequently,  results  for  stresses  are 
generally  less  accurate  than  displacements  in  any  finite 
element  solution.  This  dictates  that,  when  possible, 
extrapolation  techniques  should  be  applied  to  displace¬ 
ments  rather  than  stresses.  A  method  analogous  to  the 
stress-based  extrapolation  technique  exists  for  relative 
crack  face  displacements  [10,18].  However,  this  method 
requires  that  the  two  crack  faces  be  initially  coincident, 
so  that  this  technique  cannot  be  easily  applied  to  the 
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Table  1 


Comparison  of  numerical  modelling  and  photoelastic  results  for  stress  intensity  factors 


Applied  stress  (kPa) 

Crack  orientation  (degrees) 

Magnitudes  (kPa  m1/2) 
Computational 

Experimental 

Phase  angles  (degrees) 
Computational 

Experimental 

60.3 

0 

15.8 

19.0 

4.0 

0.0 

96.5 

15 

24.0 

30.2 

11.5 

7.8 

96.5 

30 

20.8 

20.3 

19.5 

17.1 

48.3 

45 

7.2 

8.3 

25.8 

30.0 
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current  geometry,  which  incorporates  crack  faces  with  finite 
separation. 

In  this  work,  J  integral  values  were  used  to  evaluate  K , 
and  bond  line  traction  data  was  used  to  determine  the 
phase  angle  of  K.  Results  are  presented  and  compared  with 
the  photoelasticity  results  next. 

4.  Results 

Fig.  4  shows  a  set  of  representative  contour  plots;  the 
plots  are  for  a  crack  orientation  of  45°.  The  relative  size 
of  the  two  sets  of  fringe  loops  is  caused  by  mismatch  of 
the  two  elastic  moduli;  the  shape  of  the  cryy  and  crxy  contours 
remains  invariant  with  respect  to  mode  mixity.  Other  stress 
variables,  however,  such  as  the  maximum  in-plane  stress, 
are  strong  functions  of  mode  mixity.  Table  1  summarizes 
the  results  for  both  the  magnitude  and  phase  of  K.  Experi¬ 
mental  and  computational  results  agreed  well.  With  the  T  — 
0°  and  T  =  15°  specimens,  the  photoelastic  results  were  high 
because  of  the  effect  previously  mentioned:  residual  and 
machining  stresses  incorporated  into  the  specimen  during 
fabrication  influence  the  gradient  of  the  fringe  patterns  in 
the  vertical  direction,  and  this  gradient  is  used  in  the  deter¬ 
mination  of  K. 

Table  1  also  shows  the  phase  angle  results.  The  worst 
discrepancy  between  the  numerical  and  experimental  results 
was  for  the  T  =  45°  bimaterial  specimen,  in  which  case  the 
values  differed  by  4.2°.  This  is  accurate  enough  for  most 
applications,  however,  if  necessary,  other  methods  may  be 
used  to  obtain  additional  improvements  in  accuracy  of  phase 
angle  determination  [16,17].  However,  the  method  shown 
above  is  often  preferable  because  of  its  simplicity  compared 
with  other  procedures. 

5.  Conclusions 

The  fracture  mechanics  of  interfacial  cracks  in  incom¬ 
pressible  bimaterials  subject  to  plane  strain  conditions 
closely  resembles  those  for  homogeneous  materials  subject 
to  mixed-mode  loading.  Static  indeterminacy  can  cause 
inaccuracies  unless  a  mixed  formulation  is  used.  Numerical 
modelling  can  be  used  to  characterize  the  complex  stress 
intensity  factors  for  various  specimen  geometries.  The 
magnitude  of  K  can  be  evaluated  from  J  integral  calcula¬ 
tions,  and  the  phase  angle  of  K  can  be  determined  by  extra¬ 
polating  field  variables  such  as  the  bond  line  traction  data. 
Mesh  refinements  should  include  numerous  quarter  point 
elements  near  the  crack  tip  so  that  phase  angles  can  be 
evaluated  using  extrapolation  of  bond  line  traction  data. 
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